function Data = Hydro2(Nx,Ny,tf,Lx,Ly)
% Fourier Basis in Spatial Directions
%
  close all;
  fig = figure; 
  figmetric = figure;
  load('cm', 'mycmap');
  set(fig, 'Colormap', mycmap);
  set(figmetric, 'Colormap', mycmap);
  CMAP_OFFSET = 512;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%% Initializations and Parameters                  %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PlotData = 4;
delay = 1;
dt = 1.e-6;

global Dux uxx uxy;
global Duy uyx uyy;
global DT  Tx  Ty;
global DP  Px  Py;
global DB  Bx  By;
global cc;

cc = ((sqrt(3)*pi) -(9*log(3)) +18)/18;


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%           Cash-Carp RK45 Parameters             %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C1 = 0;
C2 = 1/5;
C3 = 3/10;
C4 = 3/5;
C5 = 1;
C6 = 7/8;

A21 = 1/5;
A31 = 3/40;         A32 = 9/40;
A41 = 3/10;         A42 = -9/10;    A43 = 6/5;
A51 = -11/54;       A52 = 5/2;      A53 = -70/27;       A54 = 35/27;
A61 = 1631/55296;   A62 = 175/512;  A63 = 575/13824;    A64 = 44275/110592; A65 = 253/4096;

B1 = 37/378;        B2 = 0;         B3 = 250/621;       B4 = 125/594;       B5 = 0;         B6 = 512/1771;
b1 = 2825/27648;    b2 = 0;         b3 = 18575/48384;   b4 = 13525/55296;   b5 = 277/14336; b6 = 1/4;
E1 = B1 - b1;       E2 = B2 - b2;   E3 = B3 - b3;       E4 = B4 - b4;       E5 = B5 - b5;   E6 = B6 - b6;

aTol = 1.e-6;
eTol = 1.5d-6;
Fudge = 0.9;



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% **
    %%%%%%%%%    Build the Coordinate and Physical Grids      %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = 2*pi*(0:Nx-1)'/Nx;
y = 2*pi*(0:Ny-1)'/Ny;

[X,Y]=ndgrid(x,y);

  % Build the physical grid coordinate matrices
Xp = (Lx/2*pi)*X ;
Yp = (Ly/2*pi)*Y ;


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%         Define Helper Macro Variables           %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%                 Reserve Storeage                %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ut   = ones(Nx,Ny);
utx  = zeros(Nx,Ny);
uty  = zeros(Nx,Ny);
Dut  = zeros(Nx,Ny);

ux   = zeros(Nx,Ny);
uxx  = zeros(Nx,Ny);
uxy  = zeros(Nx,Ny);
kux  = zeros(Nx,Ny,6);
Dux  = zeros(Nx,Ny);
Eux  = zeros(Nx,Ny);

uy   = zeros(Nx,Ny);
uyx  = zeros(Nx,Ny);
uyy  = zeros(Nx,Ny);
kuy  = zeros(Nx,Ny,6);
Duy  = zeros(Nx,Ny);
Euy  = zeros(Nx,Ny);

T   = zeros(Nx,Ny);
Tx  = zeros(Nx,Ny);
Ty  = zeros(Nx,Ny);
kT  = zeros(Nx,Ny,6);
DT  = zeros(Nx,Ny);
ET  = zeros(Nx,Ny);

P   = zeros(Nx,Ny);
Px  = zeros(Nx,Ny);
Py  = zeros(Nx,Ny);
kP  = zeros(Nx,Ny,6);
DP  = zeros(Nx,Ny);
EP  = zeros(Nx,Ny);


B   = zeros(Nx,Ny);
Bx  = zeros(Nx,Ny);
By  = zeros(Nx,Ny);
kB  = zeros(Nx,Ny,6);
DB  = zeros(Nx,Ny);
EB  = zeros(Nx,Ny);


datum = zeros(Nx,Ny,5);


  global r;
  rmag = 5;
  R_SIZE = 32;
  r = rmag * tan([0:R_SIZE-1]*pi/2/R_SIZE);
  rticks = R_SIZE*2/pi*atan(linspace(r(1), r(R_SIZE), 10)/rmag);
  rticks = transpose(linspace(0,R_SIZE,6));
  rticklabels = num2str( rmag * tan(rticks*pi/2/R_SIZE), '%.g');
  gdtdt = zeros(Nx, Ny, R_SIZE);
  gdtdx = zeros(Nx, Ny, R_SIZE);
  gdtdy = zeros(Nx, Ny, R_SIZE);
  gdxdx = zeros(Nx, Ny, R_SIZE);
  gdxdy = zeros(Nx, Ny, R_SIZE);
  gdydy = zeros(Nx, Ny, R_SIZE);
  gdrdt = zeros(Nx, Ny, R_SIZE);
  gdrdx = zeros(Nx, Ny, R_SIZE);
  gdrdy = zeros(Nx, Ny, R_SIZE);




    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%      Set Variables to Initial Conditions        %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    % Set all variables to zero, padded into sufficient space


t = 0;
iteration = 1;

%T = 0*X + 5;

% for m = 3:9
%   for h = 1:10
%     phase = 2*pi*rand(1);
%     amp = 0.001*rand(1);
%     T = T + amp*cos(m*Y + phase*cos(2*X));
%     phase = 2*pi*rand(1);
%     amp = 0.001*rand(1);
%     T = T + amp*cos(m*X + phase*sin(2*Y));
%   end
% end


cx = 0.5 * Nx + 0.5;
cy = 0.5 * Ny + 0.5;
  gaussian_size = 0.1;
  sig2 = (gaussian_size*Nx).^2;
for x = 1:Nx
  for y = 1:Ny
    ux(x,y) = 0.1.*(x-cx).*exp(-((x-cx)^2 + (y-cy)^2)/sig2);
    uy(x,y) = 0.1.*(y-cy).*exp(-((x-cx)^2 + (y-cy)^2)/sig2);
  end
end
T = T + 1;



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%     Time-Evolve via Chebyshev/FFT Derivatives   %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

datum(:,:,1) =  ux(:,:);
datum(:,:,2) =  uy(:,:);
datum(:,:,3) =  T(:,:);
datum(:,:,4) =  P(:,:);
datum(:,:,5) =  B(:,:);


tic;




while t<tf
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%             Run the Cash-Carp RK45              %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    GetEOMs(ux, ...
            uy, ...
            T, ...
            P, ...
            B, ...
            t  + C1*dt,               Nx,Ny,Lx,Ly);
    kux(:,:,1) = Dux;
    kuy(:,:,1) = Duy;
    kT(:,:, 1) = DT;
    kP(:,:, 1) = DP;
    kB(:,:, 1) = DB;

    % derivatives computed here, run the metric code:
    ut = sqrt(1+ux.^2+uy.^2);
    Dut = (ux.*Dux+uy.*Duy)./ut;
    for ix = 1:Nx
    %Dy
      uty(ix,:) = (2*pi/Ly)*fourdifft(ut(ix,:),1);
    end
    for iy = 1:Ny
    %Dx
      utx(:,iy) = (2*pi/Lx)*fourdifft(ut(:,iy),1);
    end
     %[gdtdt gdtdx gdtdy gdxdx gdxdy gdydy gdrdt gdrdx gdrdy] = metric(T, ut, ux, uy, Dut, Dux, Duy, utx, uty, uxx, uxy, uyx, uyy);

    GetEOMs(ux + A21*dt*kux(:,:,1), ...
            uy + A21*dt*kuy(:,:,1), ...
            T  + A21*dt*kT(:,:,1), ...
            P  + A21*dt*kP(:,:,1), ...
            B  + A21*dt*kB(:,:,1), ...
            t  + C2*dt,               Nx,Ny,Lx,Ly);
    kux(:,:,2) = Dux;
    kuy(:,:,2) = Duy;
    kT(:,:, 2) = DT;
    kP(:,:, 2) = DP;
    kB(:,:, 2) = DB;

    GetEOMs(ux + A31*dt*kux(:,:,1) + A32*dt*kux(:,:,2), ...
            uy + A31*dt*kuy(:,:,1) + A32*dt*kuy(:,:,2), ...
            T  + A31*dt*kT(:,:,1)  + A32*dt*kT(:,:,2), ...
            P  + A31*dt*kP(:,:,1)  + A32*dt*kP(:,:,2), ...
            B  + A31*dt*kB(:,:,1)  + A32*dt*kB(:,:,2), ...
            t  + C3*dt,               Nx,Ny,Lx,Ly);
    kux(:,:,3) = Dux;
    kuy(:,:,3) = Duy;
    kT(:,:, 3) = DT;
    kP(:,:, 3) = DP;
    kB(:,:, 3) = DB;

    GetEOMs(ux + A41*dt*kux(:,:,1) + A42*dt*kux(:,:,2) + A43*dt*kux(:,:,3), ...
            uy + A41*dt*kuy(:,:,1) + A42*dt*kuy(:,:,2) + A43*dt*kuy(:,:,3), ...
            T  + A41*dt*kT(:,:,1)  + A42*dt*kT(:,:,2)  + A43*dt*kT(:,:,3), ...
            P  + A41*dt*kP(:,:,1)  + A42*dt*kP(:,:,2)  + A43*dt*kP(:,:,3), ...
            B  + A41*dt*kB(:,:,1)  + A42*dt*kB(:,:,2)  + A43*dt*kB(:,:,3), ...
            t  + C4*dt,               Nx,Ny,Lx,Ly);
    kux(:,:,4) = Dux;
    kuy(:,:,4) = Duy;
    kT(:,:, 4) = DT;
    kP(:,:, 4) = DP;
    kB(:,:, 4) = DB;

    GetEOMs(ux + A51*dt*kux(:,:,1) + A52*dt*kux(:,:,2) + A53*dt*kux(:,:,3) + A54*dt*kux(:,:,4),  ...
            uy + A51*dt*kuy(:,:,1) + A52*dt*kuy(:,:,2) + A53*dt*kuy(:,:,3) + A54*dt*kuy(:,:,4),  ...
            T  + A51*dt*kT(:,:,1)  + A52*dt*kT(:,:,2)  + A53*dt*kT(:,:,3)  + A54*dt*kT(:,:,4),   ...
            P  + A51*dt*kP(:,:,1)  + A52*dt*kP(:,:,2)  + A53*dt*kP(:,:,3)  + A54*dt*kP(:,:,4),   ...
            B  + A51*dt*kB(:,:,1)  + A52*dt*kB(:,:,2)  + A53*dt*kB(:,:,3)  + A54*dt*kB(:,:,4),   ...
            t  + C5*dt,               Nx,Ny,Lx,Ly);
    kux(:,:,5) = Dux;
    kuy(:,:,5) = Duy;
    kT(:,:, 5) = DT;
    kP(:,:, 5) = DP;
    kB(:,:, 5) = DB;

    GetEOMs(ux + A61*dt*kux(:,:,1) + A62*dt*kux(:,:,2) + A63*dt*kux(:,:,3) + A64*dt*kux(:,:,4) + A65*dt*kux(:,:,5), ...
            uy + A61*dt*kuy(:,:,1) + A62*dt*kuy(:,:,2) + A63*dt*kuy(:,:,3) + A64*dt*kuy(:,:,4) + A65*dt*kuy(:,:,5), ...
            T  + A61*dt*kT(:,:,1)  + A62*dt*kT(:,:,2)  + A63*dt*kT(:,:,3)  + A64*dt*kT(:,:,4)  + A65*dt*kT(:,:,5), ...
            P  + A61*dt*kP(:,:,1)  + A62*dt*kP(:,:,2)  + A63*dt*kP(:,:,3)  + A64*dt*kP(:,:,4)  + A65*dt*kP(:,:,5), ...
            B  + A61*dt*kB(:,:,1)  + A62*dt*kB(:,:,2)  + A63*dt*kB(:,:,3)  + A64*dt*kB(:,:,4)  + A65*dt*kB(:,:,5), ...
            t  + C6*dt,               Nx,Ny,Lx,Ly);
    kux(:,:,6) = Dux;
    kuy(:,:,6) = Duy;
    kT(:,:, 6) = DT;
    kP(:,:, 6) = DP;
    kB(:,:, 6) = DB;

    
    Eux = dt*( E1*kux(:,:,1) +E2*kux(:,:,2) +E3*kux(:,:,3) +E4*kux(:,:,4) +E5*kux(:,:,5) +E6*kux(:,:,6) );
    Euy = dt*( E1*kuy(:,:,1) +E2*kuy(:,:,2) +E3*kuy(:,:,3) +E4*kuy(:,:,4) +E5*kuy(:,:,5) +E6*kuy(:,:,6) );
    ET  = dt*( E1*kT(:,:,1)  +E2*kT(:,:,2)  +E3*kT(:,:,3)  +E4*kT(:,:,4)  +E5*kT(:,:,5)  +E6*kT(:,:,6)  );
    EP  = dt*( E1*kP(:,:,1)  +E2*kP(:,:,2)  +E3*kP(:,:,3)  +E4*kP(:,:,4)  +E5*kP(:,:,5)  +E6*kP(:,:,6)  );
    EB  = dt*( E1*kB(:,:,1)  +E2*kB(:,:,2)  +E3*kB(:,:,3)  +E4*kB(:,:,4)  +E5*kB(:,:,5)  +E6*kB(:,:,6)  );
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%             Estimate the Abs Error              %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    MaxErr = max([max(abs(Eux(:))) max(abs(Euy(:))) max(abs(ET(:))) max(abs(EP(:))) max(abs(EB(:)))]);
    
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%          If Tolerable, Update & Store           %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if MaxErr <= eTol
        ux = ux + dt*( B1*kux(:,:,1) +B2*kux(:,:,2) +B3*kux(:,:,3) +B4*kux(:,:,4) +B5*kux(:,:,5) +B6*kux(:,:,6) );
        uy = uy + dt*( B1*kuy(:,:,1) +B2*kuy(:,:,2) +B3*kuy(:,:,3) +B4*kuy(:,:,4) +B5*kuy(:,:,5) +B6*kuy(:,:,6) );
        T  = T  + dt*( B1*kT(:,:,1)  +B2*kT(:,:,2)  +B3*kT(:,:,3)  +B4*kT(:,:,4)  +B5*kT(:,:,5)  +B6*kT(:,:,6) );
        P  = P  + dt*( B1*kP(:,:,1)  +B2*kP(:,:,2)  +B3*kP(:,:,3)  +B4*kP(:,:,4)  +B5*kP(:,:,5)  +B6*kP(:,:,6) );
        B  = B  + dt*( B1*kB(:,:,1)  +B2*kB(:,:,2)  +B3*kB(:,:,3)  +B4*kB(:,:,4)  +B5*kB(:,:,5)  +B6*kB(:,:,6) );
        t  = t  +  dt;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%             Store Data As Requested             %%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


        datum(:,:,1) =  ux(:,:);
        datum(:,:,2) =  uy(:,:);
        datum(:,:,3) =  T(:,:);
        datum(:,:,4) =  P(:,:);
        datum(:,:,5) =  B(:,:);

        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%                   Plot Data                     %%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if PlotData && (mod(iteration,delay)==0)
            figure(1);
%            clf;
            surf(Xp,Yp,T);
            %image((T-1/pi)*1280 + 128);
            %zlim([1/pi-0.02 1/pi+0.02]);
            %caxis([1/pi-0.02 1/pi+0.02]);
            zlim([1-1 1+1])
            caxis([1-1 1+1])
            title(['T(t) at t = ',num2str(t)]);
            %shading interp;
    %    pmetric();
            refresh();

%             figure(2);
% %            clf;
%             ut = sqrt(1+ux.^2+uy.^2);
%             quiver(Xp,Yp,ux,uy);
%             hold on;
%             surf(Xp,Yp,ut)
%             zlim([1 1.02])
%             shading interp;
%             colormap(jet);
%             alpha(0.7);
%             title(['<O> at t = ',num2str(t)]);
%             hold off;
%             refresh();
% 
%             figure(1);
%             refresh();
        end


        % plot the metric stuff
       

        elap = toc;
        disp(['Iteration ',num2str(iteration),' at t = ',num2str(t),'s with dt = ',num2str(dt),'s took ',num2str(elap),'s']);
        iteration = iteration +1;
        ratio = Fudge * (aTol/MaxErr)^(1/5);
        if     ratio > 5
               ratio = 5;
        elseif ratio < 0.2
               ratio = 0.2;
        end
        dt = ratio * dt; 
        tic;
        
    else
        elap = toc;
        disp(['Iteration ',num2str(iteration),' FAILED with dt = ',num2str(dt),' after ',num2str(elap),'s']);
        tic;
        ratio =Fudge * (aTol/MaxErr)^(1/4);
        if ratio < 0.1
            ratio = 0.1;
        end
        dt = ratio * dt;
        disp(['Trying again with dt = ',num2str(dt)]);
    end
    
    if t + dt >= tf
        dt=tf-t;
    end
        
end



  function pmetric()
    figure(figmetric);
    plotsx = 3; 
    plotsy = 3;

    plotmag = 100;

    mdata = reshape(gdtdt(Nx/2, :,:), Ny, R_SIZE);
    subplot(plotsy, plotsx, 1);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gdtdt');

    mdata = reshape(gdtdx(Nx/2, :,:), Ny, R_SIZE);
    subplot(plotsy, plotsx, 2);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gdtdx');

    mdata = reshape(gdtdy(Nx/2, :,:), Ny, R_SIZE);
    subplot(plotsy, plotsx, 3);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gdtdy');

    mdata = reshape(gdxdx(Nx/2, :,:), Ny, R_SIZE);
    subplot(plotsy, plotsx, 4);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gdxdx');

    mdata = reshape(gdxdy(Nx/2, :,:), Ny, R_SIZE);
    subplot(plotsy, plotsx, 5);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gdxdy');

    mdata = reshape(gdydy(Nx/2, :,:), Ny, R_SIZE);
    subplot(plotsy, plotsx, 6);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gdydy');

    mdata = reshape(gdrdt(Nx/2, :,:), Ny, R_SIZE);
    subplot(plotsy, plotsx, 7);
    image(mdata*CMAP_OFFSET/plotmag+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gdrdt');

    mdata = reshape(gdrdx(Nx/2, :,:), Ny, R_SIZE);
    subplot(plotsy, plotsx, 8);
    image(mdata*CMAP_OFFSET+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gdrdx');

    mdata = reshape(gdrdy(Nx/2, :,:), Ny, R_SIZE);
    subplot(plotsy, plotsx, 9);
    image(mdata*CMAP_OFFSET+CMAP_OFFSET);
    set(gca, 'XTick', rticks);
    set(gca, 'XTickLabel', rticklabels);
    title('gdrdy');
    %if (Nx >= 64)
    %  truesize(figmetric);
    %end
  end



end




%     GetEOMs(At,Ax,Ay,pR,pI,Nx,Ny,Nz,a,b,c,Xp,Yp,IZ);
function GetEOMs(ux,uy,T,P,B,t,Nx,Ny,Lx,Ly)

% Declare Globals to Save Memory
global uxx uxy Dux;
global uyx uyy Duy;
global Tx  Ty  DT;
global Px  Py  DP;
global Bx  By  DB;
global cc;



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%   Compute Time Derivatives via Chebyshev/FFT    %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
for ix = 1:Nx
    %Dy
    uxy(ix,:) = (2*pi/Ly)*fourdifft(ux(ix,:),1);
    uyy(ix,:) = (2*pi/Ly)*fourdifft(uy(ix,:),1);
    Ty(ix,:)  = (2*pi/Ly)*fourdifft( T(ix,:),1);
    Py(ix,:)  = (2*pi/Ly)*fourdifft( P(ix,:),1);
    By(ix,:)  = (2*pi/Ly)*fourdifft( B(ix,:),1);
end
for iy = 1:Ny
    %Dx
    uxx(:,iy) = (2*pi/Lx)*fourdifft(ux(:,iy),1);
    uyx(:,iy) = (2*pi/Lx)*fourdifft(uy(:,iy),1);
    Tx(:,iy)  = (2*pi/Lx)*fourdifft( T(:,iy),1);
    Px(:,iy)  = (2*pi/Lx)*fourdifft( P(:,iy),1);
    Bx(:,iy)  = (2*pi/Lx)*fourdifft( B(:,iy),1);
end

%ut = sqrt(1+ux.^2 +uy.^2);

m11 = (2+ux.^2+uy.^2).*((-3).*T.^3.*ux+(-4).*B.*ux.*(1+uy.^2)+(-2).*P.*uy.*(2+(-1).*ux.^2+uy.^2));
m12 = (-1).*(2+ux.^2+uy.^2).*((3.*T.^3+(-4).*B.*(1+ux.^2)).*uy+2.*P.*ux.*(2+ux.^2+(-1).*uy.^2));
m13 = (-6).*T.^2.*(1+ux.^2+uy.^2).*(2+ux.^2+uy.^2);
%m14 = 0;
%m15 = 0;

m21 = (1/2).*(1+ux.^2+uy.^2).^(-1).*(2+ux.^2+uy.^2).*(3.*T.^3.*(1+ux.^2+ ...
  uy.^2).*(2+ux.^2+2.*uy.^2)+2.*P.*ux.*uy.*(2+3.*(ux.^2+ux.^4)+ ...
  uy.^2+(-3).*uy.^4)+4.*B.*(ux.^4.*(1+(-3).*uy.^2)+2.*(1+uy.^2)+(-3) ...
  .*ux.^2.*((-1)+uy.^4)));
m22 = (-1/2).*(1+ux.^2+uy.^2).^(-1).*(2+ux.^2+uy.^2).*(ux.*uy.*(3.* ...
  T.^3.*(1+ux.^2+uy.^2)+(-4).*B.*(1+ux.^2).*(1+3.*ux.^2+3.*uy.^2))+ ...
  2.*P.*((-4)+3.*ux.^4+3.*ux.^6+(-6).*uy.^2+(-4).*uy.^4+(-1).* ...
  ux.^2.*(4+3.*(uy.^2+uy.^4))));
%m23 = 0;
m24 = 2.*uy.*(2+ux.^2+uy.^2).*(2+3.*ux.^2+uy.^2);
m25 = 4.*(ux+ux.^3).*(2+ux.^2+uy.^2);

m31 = (1/2).*(1+ux.^2+uy.^2).^(-1).*(2+ux.^2+uy.^2).*((-1).*ux.*uy.*(3.* ...
  T.^3.*(1+ux.^2+uy.^2)+4.*B.*(1+uy.^2).*(1+3.*ux.^2+3.*uy.^2))+2.* ...
  P.*(3.*ux.^2.*(2+uy.^2)+ux.^4.*(4+3.*uy.^2)+(-1).*(1+uy.^2).*((-4) ...
  +3.*uy.^4)));
m32 = (1/2).*(1+ux.^2+uy.^2).^(-1).*(2+ux.^2+uy.^2).*(3.*T.^3.*(1+ux.^2+ ...
  uy.^2).*(2+2.*ux.^2+uy.^2)+4.*B.*((-2).*(1+ux.^2)+3.*((-1)+ux.^4) ...
  .*uy.^2+((-1)+3.*ux.^2).*uy.^4)+2.*P.*ux.*uy.*(2+ux.^2+(-3).* ...
  ux.^4+3.*(uy.^2+uy.^4)));
%m33 = 0;
m34 = 2.*ux.*(2+ux.^2+uy.^2).*(2+ux.^2+3.*uy.^2);
m35 = (-4).*(2+ux.^2+uy.^2).*(uy+uy.^3);


m41 = (-1/3).*T.^(-4).*(2+ux.^2+uy.^2).*(3.*T.^3.*uy.*(2.*P.*ux.*(ux+( ...
  -1).*uy).*uy.*(ux+uy)+T.^3.*(1+uy.^2))+4.*B.*(1+uy.^2).*((-3).* ...
  T.^3.*(1+ux.^2).*uy+cc.*ux.*((-3).*T.^3.*ux.*uy+2.*P.*(2+ux.^2+ ...
  uy.^2)))+cc.*P.*(4.*P.*uy.*((-1).*ux.^4+(2+uy.^2).^2)+3.*T.^3.* ...
  ux.*(6+uy.^2+(-2).*uy.^4+ux.^2.*(5+2.*uy.^2))));
m42 = (-1/3).*T.^(-4).*(2+ux.^2+uy.^2).*(3.*T.^3.*ux.*(T.^3.*(1+ux.^2)+ ...
  4.*B.*(1+ux.^2).*(1+uy.^2)+2.*P.*ux.*uy.*((-1).*ux.^2+uy.^2))+cc.* ...
  (12.*B.*T.^3.*ux.*(1+ux.^2).*uy.^2+4.*P.^2.*ux.*((2+ux.^2).^2+(-1) ...
  .*uy.^4)+P.*uy.*((-8).*B.*(1+ux.^2).*(2+ux.^2+uy.^2)+3.*T.^3.*(6+ ...
  ux.^2+(-2).*ux.^4+(5+2.*ux.^2).*uy.^2))));
%m43 = 0;
m44 = (-2).*cc.*T.^(-1).*(1+ux.^2+uy.^2).*(2+ux.^2+uy.^2).^2;
%m45 = 0;


m51 = (1/3).*T.^(-4).*(2+ux.^2+uy.^2).*(16.*B.^2.*cc.*ux.*(1+uy.^2).*(2+ ...
  ux.^2+uy.^2)+3.*T.^3.*(T.^3.*ux.*(2+ux.^2+3.*uy.^2)+2.*P.*uy.*((1+ ...
  cc).*ux.^4+(2+uy.^2).*(2+(1+cc).*uy.^2)+(-2).*ux.^2.*((-2)+cc+(1+ ...
  cc).*uy.^2)))+2.*B.*(6.*T.^3.*ux.*uy.^2.*((-1).*ux.^2+uy.^2)+cc.*( ...
  3.*T.^3.*ux.*(6+7.*uy.^2+2.*uy.^4+ux.^2.*(3+(-2).*uy.^2))+4.*P.* ...
  uy.*((-1).*ux.^4+(2+uy.^2).^2))));
m52 = (-1/3).*T.^(-4).*(2+ux.^2+uy.^2).*(16.*B.^2.*cc.*(1+ux.^2).*uy.*( ...
  2+ux.^2+uy.^2)+2.*B.*((-4).*cc.*P.*ux.*(2+ux.^2).^2+(-3).*T.^3.*( ...
  6.*cc+7.*cc.*ux.^2+2.*(1+cc).*ux.^4).*uy+3.*T.^3.*((-3).*cc+2.*(1+ ...
  cc).*ux.^2).*uy.^3+4.*cc.*P.*ux.*uy.^4)+3.*T.^3.*(T.^3.*uy.*(2+3.* ...
  ux.^2+uy.^2)+2.*P.*ux.*((2+ux.^2).*(2+(1+cc).*ux.^2)+(-2).*((-2)+ ...
  cc+(1+cc).*ux.^2).*uy.^2+(1+cc).*uy.^4)));
%m53 = 0;
%m54 = 0;
m54 = 4.*cc.*T.^(-1).*(1+ux.^2+uy.^2).*(2+ux.^2+uy.^2).^2;






b1 = (1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).*(6.*T.^2.*(Tx.*ux+Ty.* ...
  uy).*(1+ux.^2+uy.^2)+3.*T.^3.*(1+ux.^2+uy.^2).*(uxx+uyy)+2.*(2.* ...
  P.*uxy+2.*B.*((1+ux.^2).*uxx.*(1+uy.^2)+ux.*uxy.*(uy+uy.^3)+(-1).* ...
  (ux+ux.^3).*uy.*uyx+(-1).*(1+ux.^2).*(1+uy.^2).*uyy)+P.*(2.*uyx+ ...
  ux.^4.*uyx+uy.^2.*(uxy.*(3+uy.^2)+uyx)+ux.^2.*(uxy+(-1).*uxy.* ...
  uy.^2+3.*uyx+(-1).*uy.^2.*uyx)+ux.*uy.^3.*(uxx+(-1).*uyy)+ux.^3.* ...
  uy.*((-1).*uxx+uyy))));
b2 = (-1/2).*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).*(6.*T.^2.*Tx.*( ...
  1+ux.^2+uy.^2)+8.*Bx.*(1+ux.^2).*(1+ux.^2+uy.^2)+8.*Px.*ux.*uy.*( ...
  1+ux.^2+uy.^2)+4.*Py.*(1+ux.^2+uy.^2).*(2+ux.^2+uy.^2)+uxx.*(3.* ...
  T.^3.*ux.*(1+ux.^2+uy.^2)+4.*B.*ux.*(1+uy.^2+ux.^2.*(1+(-3).* ...
uy.^2))+2.*P.*uy.*(ux.^2.*(4+3.*ux.^2+(-3).*uy.^2)+4.*(1+uy.^2)))+ ...
  2.*uxy.*(P.*ux.*((-2)+ux.^2+((-5)+3.*ux.^2).*uy.^2+(-3).*uy.^4)+ ...
  3.*uy.*((-2).*B.*ux.^2.*(1+uy.^2)+T.^3.*(1+ux.^2+uy.^2)))+(-2).* ...
  ux.*((-6).*B.*ux.*(1+ux.^2).*uy+P.*(2+3.*ux.^4+(-1).*uy.^2+ux.^2.* ...
  (5+(-3).*uy.^2))).*uyx+(12.*B.*(ux+ux.^3).*(1+uy.^2)+(-3).*T.^3.* ...
  ux.*(1+ux.^2+uy.^2)+2.*P.*uy.*(4+4.*ux.^2+(-3).*ux.^4+(4+3.*ux.^2) ...
  .*uy.^2)).*uyy);
b3 = (-1/2).*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).*(6.*T.^2.*Ty.*( ...
  1+ux.^2+uy.^2)+8.*Py.*ux.*uy.*(1+ux.^2+uy.^2)+(-8).*By.*(1+uy.^2) ...
  .*(1+ux.^2+uy.^2)+4.*Px.*(1+ux.^2+uy.^2).*(2+ux.^2+uy.^2)+(-2).* ...
  uxy.*uy.*(6.*B.*ux.*uy.*(1+uy.^2)+P.*(2+5.*uy.^2+3.*uy.^4+(-1).* ...
  ux.^2.*(1+3.*uy.^2)))+uxx.*((-3).*uy.*(4.*B.*(1+ux.^2).*(1+uy.^2)+ ...
  T.^3.*(1+ux.^2+uy.^2))+2.*P.*ux.*(4+4.*uy.^2+(-3).*uy.^4+ux.^2.*( ...
  4+3.*uy.^2)))+2.*(3.*T.^3.*ux.*(1+ux.^2+uy.^2)+uy.*(6.*B.*ux.*(1+ ...
  ux.^2).*uy+P.*((-2)+uy.^2+ux.^2.*((-5)+(-3).*ux.^2+3.*uy.^2)))).* ...
  uyx+(-1).*((-3).*T.^3.*uy.*(1+ux.^2+uy.^2)+4.*B.*uy.*(1+ux.^2+ ...
  uy.^2+(-3).*ux.^2.*uy.^2)+2.*P.*ux.*((-4)+(-4).*uy.^2+(-3).*uy.^4+ ...
  ux.^2.*((-4)+3.*uy.^2))).*uyy);
b4 = (-1/3).*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).*((-12).*cc.*Px.* ...
  T.^(-1).*ux+(-1).*T.^(-4).*(4.*cc.*P.^2.*(2+ux.^2+uy.^2).*(uxy.*( ...
  2+ux.^2+(-1).*((-3)+ux.^2).*uy.^2+uy.^4)+ux.^4.*uyx+(-1).*ux.^2.*( ...
  (-3)+uy.^2).*uyx+(2+uy.^2).*uyx+ux.*uy.^3.*(uxx+(-1).*uyy)+ux.^3.* ...
  uy.*((-1).*uxx+uyy))+3.*T.^3.*(T.^3.*(1+ux.^2+uy.^2).*(uxy+uxy.* ...
  uy.^2+uyx+ux.^2.*uyx)+4.*B.*(1+ux.^2).*(1+uy.^2).*((-1).*uxy.*(1+ ...
  uy.^2)+uyx+ux.^2.*uyx+ux.*uy.*((-1).*uxx+uyy))+2.*cc.*(Py.*uy.*(1+ ...
  ux.^2+uy.^2).*(2+ux.^2+uy.^2)+Px.*ux.*(ux.^2+uy.^2).*(3+ux.^2+ ...
  uy.^2)+2.*B.*((-1).*(1+uy.^2).*(ux.^3.*uxx.*uy+uxy.*((-1)+ux.^2.*( ...
  (-1)+uy.^2)))+(1+ux.^2).*((-1)+((-1)+ux.^2).*uy.^2).*uyx+ux.*(1+ ...
  ux.^2).*uy.^3.*uyy)))+P.*(6.*T.^4.*(1+ux.^2+uy.^2).*(2+ux.^2+ ...
  uy.^2)+(-8).*B.*cc.*(2+ux.^2+uy.^2).*((-1).*(1+ux.^2).*uxx.*(1+ ...
  uy.^2)+(-1).*ux.*uxy.*(uy+uy.^3)+(ux+ux.^3).*uy.*uyx+(1+ux.^2).*( ...
  1+uy.^2).*uyy)+3.*T.^3.*(cc.*(uxx.*(6+9.*ux.^2+5.*ux.^4+(9+4.* ...
  ux.^2+2.*ux.^4).*uy.^2+(3+(-2).*ux.^2).*uy.^4)+(-2).*ux.*(ux+(-1) ...
  .*uy).*uy.*(ux+uy).*((-1).*uxy.*uy.^2+ux.^2.*uyx)+(6+9.*uy.^2+5.* ...
  uy.^4+ux.^4.*(3+(-2).*uy.^2)+ux.^2.*(9+4.*uy.^2+2.*uy.^4)).*uyy)+( ...
  -2).*ux.*(ux+(-1).*uy).*uy.*(ux+uy).*((-1).*uxy.*(1+uy.^2)+uyx+ ...
  ux.*((-1).*uxx.*uy+ux.*uyx+uy.*uyy))))));
b5 = (-1/3).*T.^(-4).*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+uy.^2).*((-16) ...
  .*B.^2.*cc.*(2+ux.^2+uy.^2).*((-1).*(1+ux.^2).*uxx.*(1+uy.^2)+(-1) ...
  .*ux.*uxy.*(uy+uy.^3)+(ux+ux.^3).*uy.*uyx+(1+ux.^2).*(1+uy.^2).* ...
  uyy)+2.*B.*(6.*T.^4.*(1+ux.^2+uy.^2).*(2+ux.^2+uy.^2)+4.*cc.*P.*( ...
  2+ux.^2+uy.^2).*(uxy.*(2+ux.^2+(-1).*((-3)+ux.^2).*uy.^2+uy.^4)+ ...
  ux.^4.*uyx+(-1).*ux.^2.*((-3)+uy.^2).*uyx+(2+uy.^2).*uyx+ux.* ...
  uy.^3.*(uxx+(-1).*uyy)+ux.^3.*uy.*((-1).*uxx+uyy))+3.*T.^3.*(cc.*( ...
  uxx.*(6+9.*uy.^2+3.*uy.^4+ux.^4.*(3+(-2).*uy.^2)+ux.^2.*(9+2.* ...
  uy.^2.*(5+uy.^2)))+2.*ux.*uy.*(ux.^4.*uyx+(-1).*ux.^2.*((-1)+ ...
  uy.^2).*(uxy+uyx)+uy.^2.*(uxy+uxy.*uy.^2+uyx))+(3.*(2+3.*ux.^2+ ...
  ux.^4)+(9+2.*ux.^2.*(5+ux.^2)).*uy.^2+(3+(-2).*ux.^2).*uy.^4).* ...
  uyy)+2.*ux.*(ux+(-1).*uy).*uy.*(ux+uy).*((-1).*uxy.*(1+uy.^2)+uyx+ ...
  ux.*((-1).*uxx.*uy+ux.*uyx+uy.*uyy))))+3.*T.^3.*((-8).*((-1)+cc).* ...
  P.*uxy+4.*Bx.*cc.*ux.*(1+ux.^2+uy.^2).*(2+ux.^2+uy.^2)+T.^3.*(1+ ...
  ux.^2+uy.^2).*(uxx.*(2+ux.^2+uy.^2)+2.*ux.*uy.*(uxy+(-1).*uyx)+( ...
  -1).*(2+ux.^2+uy.^2).*uyy)+2.*(2.*By.*cc.*uy.*(1+ux.^2+uy.^2).*(2+ ...
  ux.^2+uy.^2)+P.*(uxy.*uy.^2.*(8+(-4).*cc+(5+cc).*uy.^2+(1+cc).* ...
  uy.^4)+(-1).*(1+cc).*ux.^6.*uyx+((-1)+cc).*(2+uy.^2).^2.*uyx+ ...
  ux.^4.*(uxy.*(1+(-1).*cc+(1+cc).*uy.^2)+((-5)+(-1).*cc+2.*(1+cc).* ...
  uy.^2).*uyx)+ux.^2.*((-2).*uxy.*((-2)+2.*cc+(-1).*uy.^2+uy.^4+cc.* ...
  uy.^4)+(-1).*(8+2.*uy.^2+uy.^4+cc.*((-4)+uy.^4)).*uyx)+(1+cc).* ...
  ux.^5.*uy.*(uxx+(-1).*uyy)+2.*ux.^3.*uy.*((-1).*uxx.*((-2)+cc+(1+ ...
  cc).*uy.^2)+((-2)+(-1).*cc+(1+cc).*uy.^2).*uyy)+ux.*uy.*(uxx.*(2+ ...
  uy.^2).*(2+(1+cc).*uy.^2)+(-1).*(4+(-2).*((-2)+cc).*uy.^2+(1+cc).* ...
  uy.^4).*uyy)))));



Dux = (b5.*(m25.*m34.*m42+(-1).*m24.*m35.*m42+(-1).*m25.*m32.*m44+m22.* ...
  m35.*m44)+b3.*m25.*m44.*m52+(-1).*b2.*m35.*m44.*m52+b3.*m24.*m42.* ...
  m54+(-1).*b2.*m34.*m42.*m54+(-1).*b3.*m22.*m44.*m54+b2.*m32.*m44.* ...
  m54+b4.*((-1).*m25.*m34.*m52+m24.*m35.*m52+(-1).*m24.*m32.*m54+ ...
  m22.*m34.*m54)).*(m22.*m35.*m44.*m51+(-1).*m21.*m35.*m44.*m52+ ...
 m25.*(m34.*m42.*m51+(-1).*m32.*m44.*m51+(-1).*m34.*m41.*m52+m31.* ...
  m44.*m52)+m22.*m34.*m41.*m54+(-1).*m21.*m34.*m42.*m54+(-1).*m22.* ...
  m31.*m44.*m54+m21.*m32.*m44.*m54+m24.*((-1).*m35.*m42.*m51+m35.* ...
  m41.*m52+(-1).*m32.*m41.*m54+m31.*m42.*m54)).^(-1);
Duy = (b5.*(m25.*m34.*m41+(-1).*m24.*m35.*m41+(-1).*m25.*m31.*m44+m21.* ...
  m35.*m44)+((-1).*b4.*m25.*m34+b4.*m24.*m35+b3.*m25.*m44+(-1).*b2.* ...
  m35.*m44).*m51+((-1).*b4.*m24.*m31+b4.*m21.*m34+b3.*m24.*m41+(-1) ...
  .*b2.*m34.*m41+(-1).*b3.*m21.*m44+b2.*m31.*m44).*m54).*(m35.*( ...
 m24.*m42.*m51+(-1).*m22.*m44.*m51+(-1).*m24.*m41.*m52+m21.*m44.* ...
  m52)+m25.*((-1).*m34.*m42.*m51+m32.*m44.*m51+m34.*m41.*m52+(-1).* ...
  m31.*m44.*m52)+(m24.*m32.*m41+(-1).*m22.*m34.*m41+(-1).*m24.*m31.* ...
  m42+m21.*m34.*m42+m22.*m31.*m44+(-1).*m21.*m32.*m44).*m54).^(-1);
DT  = m13.^(-1).*(m35.*((-1).*m24.*m42.*m51+m22.*m44.*m51+m24.*m41.*m52+ ...
  (-1).*m21.*m44.*m52)+m25.*(m34.*m42.*m51+(-1).*m32.*m44.*m51+(-1) ...
  .*m34.*m41.*m52+m31.*m44.*m52)+((-1).*m24.*m32.*m41+m22.*m34.*m41+ ...
  m24.*m31.*m42+(-1).*m21.*m34.*m42+(-1).*m22.*m31.*m44+m21.*m32.* ...
  m44).*m54).^(-1).*(b5.*(m25.*m34+(-1).*m24.*m35).*(m12.*m41+(-1).* ...
  m11.*m42)+b5.*((-1).*m12.*m25.*m31+m11.*m25.*m32+m12.*m21.*m35+( ...
  -1).*m11.*m22.*m35).*m44+b1.*m25.*m34.*m42.*m51+(-1).*b1.*m24.* ...
  m35.*m42.*m51+b3.*m12.*m25.*m44.*m51+(-1).*b1.*m25.*m32.*m44.*m51+ ...
  (-1).*b2.*m12.*m35.*m44.*m51+b1.*m22.*m35.*m44.*m51+(-1).*b1.* ...
 m25.*m34.*m41.*m52+b1.*m24.*m35.*m41.*m52+(-1).*b3.*m11.*m25.* ...
m44.*m52+b1.*m25.*m31.*m44.*m52+b2.*m11.*m35.*m44.*m52+(-1).*b1.* ...
  m21.*m35.*m44.*m52+(-1).*b4.*(m25.*m34+(-1).*m24.*m35).*(m12.*m51+ ...
  (-1).*m11.*m52)+b4.*((-1).*m12.*m24.*m31+m11.*m24.*m32+m12.*m21.* ...
  m34+(-1).*m11.*m22.*m34).*m54+(b3.*(m12.*m24.*m41+(-1).*m11.*m24.* ...
  m42+(-1).*m12.*m21.*m44+m11.*m22.*m44)+b2.*((-1).*m12.*m34.*m41+ ...
  m11.*m34.*m42+m12.*m31.*m44+(-1).*m11.*m32.*m44)+b1.*((-1).*m24.* ...
  m32.*m41+m22.*m34.*m41+m24.*m31.*m42+(-1).*m21.*m34.*m42+(-1).* ...
  m22.*m31.*m44+m21.*m32.*m44)).*m54);
DP  = (m35.*(m24.*m42.*m51+(-1).*m22.*m44.*m51+(-1).*m24.*m41.*m52+m21.* ...
  m44.*m52)+m25.*((-1).*m34.*m42.*m51+m32.*m44.*m51+m34.*m41.*m52+( ...
  -1).*m31.*m44.*m52)+(m24.*m32.*m41+(-1).*m22.*m34.*m41+(-1).*m24.* ...
  m31.*m42+m21.*m34.*m42+m22.*m31.*m44+(-1).*m21.*m32.*m44).*m54).^( ...
  -1).*(b5.*((-1).*m25.*m32.*m41+m22.*m35.*m41+m25.*m31.*m42+(-1).* ...
  m21.*m35.*m42)+(-1).*(b3.*m25+(-1).*b2.*m35).*(m42.*m51+(-1).* ...
 m41.*m52)+((-1).*b3.*m22.*m41+b2.*m32.*m41+b3.*m21.*m42+(-1).*b2.* ...
  m31.*m42).*m54+b4.*(m25.*m32.*m51+(-1).*m22.*m35.*m51+(-1).*m25.* ...
  m31.*m52+m21.*m35.*m52+m22.*m31.*m54+(-1).*m21.*m32.*m54));
DB  = (b5.*(m24.*m32.*m41+(-1).*m22.*m34.*m41+(-1).*m24.*m31.*m42+m21.* ...
  m34.*m42+m22.*m31.*m44+(-1).*m21.*m32.*m44)+((-1).*b4.*m24.*m32+ ...
  b4.*m22.*m34+b3.*m24.*m42+(-1).*b2.*m34.*m42+(-1).*b3.*m22.*m44+ ...
  b2.*m32.*m44).*m51+(b4.*m24.*m31+(-1).*b4.*m21.*m34+(-1).*b3.* ...
 m24.*m41+b2.*m34.*m41+b3.*m21.*m44+(-1).*b2.*m31.*m44).*m52).*( ...
 m35.*(m24.*m42.*m51+(-1).*m22.*m44.*m51+(-1).*m24.*m41.*m52+m21.* ...
  m44.*m52)+m25.*((-1).*m34.*m42.*m51+m32.*m44.*m51+m34.*m41.*m52+( ...
  -1).*m31.*m44.*m52)+(m24.*m32.*m41+(-1).*m22.*m34.*m41+(-1).*m24.* ...
  m31.*m42+m21.*m34.*m42+m22.*m31.*m44+(-1).*m21.*m32.*m44).*m54).^( ...
  -1);


end





